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GUARANTEED AND SHARP A POSTERIORI ERROR ESTIMATES 
IN ISOGEOMETRIC ANALYSIS 

STEFAN K. KLEISS AND SATYENDRA K. TOMAR 

Abstract. We present functional- type a posteriori error estimates in isogeometric analysis 
(IGA). These estimates, derived on functional grounds, provide guaranteed and sharp upper 
bounds of the exact error in the energy norm. By exploiting the properties of non-uniform 
rational B-splines (NURBS), we present efficient computation of these error estimates. The 
numerical realization and the quality of the computed error distribution are addressed. The po- 
tential and the limitations of the proposed approach are illustrated using several computational 
examples. 



1. Introduction 



The geometry representations in finite element methods (FEM) and computer aided design 
(CAD) have been developed independent of each other, and are optimized for the purposes 
(S| ■ within their respective fields. As a consequence, the representations are different from each 

other, and a transfer of geometry information from CAD to FEM programmes (and vice versa) 
requires a transformation of geometry data. These transformations are, in general, not only 
costly, but also prone to approximation errors, and may require manual input. 

Isogeometric analysis (IGA), introduced by Hughes et al. [ ], see also [11], aims at closing 
this gap between FEM and CAD. The key observation is that it is a widespread standard in 
CAD to use geometry representations based on non-uniform rational B-splines (NURBS), and 
that these NURBS basis functions have properties which make them suitable as basis functions 
for FEM. Instead of transforming the geometry data to a conventional FEM representation, the 
£\| ' original geometry description is used directly, and the underlying NURBS functions are used 

as basis for the discrete solution. This way, the geometry is represented exactly in the sense 
that the geometry obtained from CAD is not changed. Thus, the need for data transformation 
r»^ ■ is eliminated, and furthermore, the exact representation from the coarsest mesh is preserved 

t^ . throughout the refinement process. IGA has been thoroughly studied and analyzed (see, e.g., 

t^ \ [3, 7, 12, 23, 36]), and its potential has been shown by successful applications to a wide range 

O ■ of problems (see, e.g., [5, 6, 10, 17, 27]). 

As mentioned above, the most widely used geometry representations in CAD are based on 
NURBS. The straightforward definition of NURBS basis functions leads to a tensor-product 
structure of the basis functions, and thus of the discretization. Since naive mesh refinement 



S^ in a tensor-product-setting has global effects, the development of local refinement strategies for 

isogeometric analysis is a subject of current active research. Such local refinement techniques 
include, for example, T-splines [4, 26, 33, 34, 35], truncated hierarchical B-splines (THB-splines) 
[21], polynomial splines over hierarchical T-meshes (PHT-splines) [ ], and locally-refineable 
splines (LR-splines) [14]. 

The issue of adaptive, local refinement is closely linked to the question of efficient a posteriori 
error estimation (see, e.g., [1, 32] for a general overview on error estimators). In the light of 
adaptive refinement, an error estimator has to identify the areas where further refinement is 
needed due to the local error being significantly larger than in the rest of the domain. Hence, an 
accurate indication of the error distribution is essential. However, a posteriori error estimation 
in isogeometric analysis is still in an infancy stage. To the best of the authors knowledge, the 
only published results are [15, 37, 38]. In [15, 37], the authors used the a posteriori error esti- 
mates based on hierarchical bases [2]. Its reliability and efficiency is subjected to the saturation 
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assumption on the (enlarged) underlying space and the constants in the strengthened Cauchy 
inequality. As the authors remarked, the first assumption is critical and its validity depends 
on the considered example. Moreover, an accurate estimation of constants in the strengthened 
Cauchy inequality requires the solution of generalized minimum eigenvalue problem. In [38], 
the authors used the residual-based a posteriori error estimates, which require the computation 
of constants in Clement-type interpolation operators. Such constants are mesh (element) de- 
pendent, often unknown or incomputable for general element shape; and the global constant 
often over-estimates the local constants, and thus the exact error. Furthermore, Zienkiewicz- 
Zhu type a posteriori error estimates are based on post-processing of approximate solutions, 
and depends on the super convergence properties of the underlying basis. To the best of authors' 
knowledge, super convergence properties for B-splines (NURBS) functions are not known. Sum- 
marily, in most general situations, the reliability and efficiency of these methods often depend 
on undetermined constants. 

In this paper, we present a functional-type a posteriori error estimator for isogeometric dis- 
cretizations. These error estimates, which were introduced in [29, 30, 31] and have been studied 
for various fields (see [32] and the references therein), provide guaranteed, fully computable and 
sharp bounds (without any generic undetermined constants). These estimates are derived on 
purely functional grounds (based on integral identities or functional analysis) and are thus ap- 
plicable to any conforming approximation in the respective space. For elliptic problems with the 
weak solution u E Hq(Q), these error bounds involve computing a free function y E i?(f2, div). 

Two aspects motivate the application of functional-type error estimates in IGA: Firstly, un- 
like the standard Lagrange basis functions, NURBS basis functions of degree p are, in general, 
globally C p-1 -continuous. Hence, NURBS basis functions of degree p > 2 are, in general, at 
least C 1 -continuous, and therefore, automatically in i7(fi,div). Thereby, we avoid construct- 
ing complicated functions in iJ(f2, div), in particular for higher degrees (see, e.g., [8, 9, 18]). 
Secondly, since the considered problem is solved in an isogeometric setting, an efficient imple- 
mentation of NURBS basis functions is readily available, which can be used to construct the 
above mentioned function y. Hence, applying the technique of functional- type a posteriori error 
estimation in a setting that relies only on the use of already available NURBS basis functions is 
greatly appealing. 

The remainder of this paper is organized as follows. In Section 2, we define the model problem, 
and recall the definition and some important properties of B-spline and NURBS basis functions. 
After briefly recalling functional-type a posteriori error estimator for IGA in Section 3, we 
propose our quality and local error indicator. In Section 4, we discuss a cost-efficient realization 
of the proposed error estimator using an illustrative numerical example. Further numerical 
examples are presented in Section 5, and finally, conclusions are drawn in Section 6. 

2. Preliminaries 

In order to fix notation and to provide an overview, we define the model problem and recall 
the definition and some aspects of isogeometric analysis in this section. 

2.1. Model Problem. Let £1 C R 2 be an open, bounded and connected Lipschitz domain with 
boundary dfl. We shall consider the following model problem: 
Find the scalar function u : Vt — >► R such that 

m -dw(AVu) = f in £}, 

^ ' u — ud on r^ = dQ, 

where A, / and gp are given data. We assume that A is a symmetric positive definite matrix 
and has a positive inverse A -1 , and that there exist constants ci, C2 > such that 

(2) ci|£| 2 <j4M<C2|e| 2 , V£GM 2 . 

Then, the norms 



(3) \\v\\\ = / Av-vdx, \\v\\ 2 A = / A x v 



• v dx, 



are equivalent to the L 2 -norm \\v\\ 2 = J^ v • v dx. The weak form of problem (1) can be written 
as follows: 

Find u <EV g , such that 

(4) a(u,v) = f(v), VveVo, 

where Vq C H 1 ^) contains the functions which vanish on T^, and V g C H 1 ^) contains the 
functions satisfying the Dirichlet boundary conditions u — up on Tjj. We assume that the 
problem data A, / and gjj are given such that the bilinear form a(-, •) is bounded, symmetric 
and positive definite, and that /(•) is a bounded linear functional. The energy norm of a function 
v is given by || Vi>||a = \Ja{y, v). Note that we have considered the Dirichlet problem only for the 
sake of simplicity. Functional-type error estimates can be easily generalized to mixed problems, 
see, e.g., [25, 32]. 

We discretize the problem (4) in the standard way by choosing a finite-dimensional space 
Vh C V g and looking for a discrete solution % G l^. This leads to a linear system of equations 
of the form 

(5) K h u h =l h , 

where K_ h is the stiffness matrix induced by the bilinear form a(-, •), / is the load vector, and 
u h is the coefficient vector of the discrete solution u^. 

2.2. B-Splines, NURBS and Isogeometric Analysis. We briefly recall the definition of 
B-spline basis functions and NURBS mappings. We only provide the basic definitions and prop- 
erties relevant for the scope of this paper. For detailed discussions of NURBS basis functions, 
geometry mappings and their properties, we refer to, e.g., [11, 12, 22, 28] and the references 
therein. The following standard definitions and statements can also be found there. 

Let p be a non-negative degree and let s = (si, . . . , s m ) be a knot vector with s^ < Si+i for 
all i. We consider only open knot vectors, i.e., knot vectors s where the multiplicity of a knot 
is at most p, except for the first and last knot which have multiplicity p + 1. For simplicity, 
we assume that si = and s m = 1, which can be easily achieved by a suitable scaling. The 
n — m — p — 1 univariate B-spline basis functions Bf : (0, 1) — > M, i = 1, . . . , n, are defined 
recursively as follows: 

1 for si < £ < s i+ i 



B *>°® ~ { else 



Si-\-p $i Si-\-p-\-l Si-\-l 

Whenever a zero denominator appears in the definition above, the corresponding function Bf 
is zero, and the whole term is considered to be zero. For open knot vectors, the first and last 
basis function are interpolatory at the first and the last knot, respectively. The derivatives of 
B-spline basis functions are given by the following formula: 

<W»(0 = T^— B t P -i(0-- ^— £f+i, P -i(0- 

^i-\-p ^i ^i-\-p-\-l &i-\-l 

B-spline basis functions of degree p are, in general, globally C p-1 -continuous. In the presence 
of multiple knots, the continuity reduces according to the multiplicity, i.e., if a knot appears k 
times, the continuity of a B-spline basis function of degree p at that knot is C p ~ k . 

Let {Bf Yl^ and {^j 9 }?=i be two families of B-spline basis functions defined by the degrees 
p and q, and the open knot vectors 

s = (si, . . . , s ni+p _|_i), t = (£i, . . . , £ n2+(? _j_i), 

respectively. We denote the set of all double-indices (i, j) by 

Zr = {(hj) ' i € {!,... ,ni},j G {l,...,n 2 }}. 



Let wuj), (hj) £ Zr, be positive weights. The bivariate NURBS basis functions i?^j) (^1,^2)5 
(i,j) E 2# are defined as follows: 

R {P P s w _m *k(ft) *yfe) 

(M) ^^ 2J E (M)e z^(M)^ )P (6)^ g fe)- 

The continuity of the B-spline basis functions is inherited by the NURBS basis functions. Note 
that B-splines can be seen as a special case of NURBS with all weights being equal to one. 
Hence, we will not distinguish between these two and we will only use the term NURBS in the 
remainder of the paper. 
The set of functions 

V h = span{ifyj), (ij) G 1r}, 

associated with the parameter domain ft = (0, l) 2 , is uniquely determined by the degrees p and 
q : the knot vectors s and t, and the weights w. Given the set of functions Vh and a control net 
of control points Puj) G K 2 , where (i, j) G 2#, the two-dimensional NURBS-surface G : Ct ^ Q 
is defined by 

(6) G(£i,&)= ^ %;)(£i,6)^)- 

We refer to fi = G(fi) as the physical domain. We assume that the geometry mapping is 
continuous and bijective (i.e., not self-penetrating), which are natural assumptions for CAD- 
applications. 

In isogeometric analysis, the isoparametric principle is applied by using the same basis func- 
tions for representing both the geometry and the discrete solution Uh- For detailed discussion, 
we refer the reader to, e.g., [11, 12, 22]. The discrete solution Uh on the physical domain Vt is 
represented as follows: 

(7) u h (x) = J2 HhJ) ( R (iJ) ° G_1 )( x )' 

(iJ)ex R 

where uuj\ G R are real- valued coefficients which form the coefficient vector u h . The discrete 
functions space is thus defined by 

V h = span{ifyj) o G" 1 , (ij) G Xr\. 

The initial mesh, and thereby the basis functions on this initial mesh, are assumed to be given 
via the geometry representation of the computational domain, i.e., the initial discretization is 
already determined by the problem domain. The exact representation of the geometry on the 
initial (coarsest) level is preserved in the process of mesh refinement. 

As mentioned in the introduction, the straightforward definition of NURBS basis functions, 
leads to a tensor-product structure of the discretization, which is the focus of this paper. Nev- 
ertheless, the error estimator presented herein is also applicable to local refinement techniques 
(e.g., T-splines, THB-splines, PHT-splines, LR-splines, see Section 1) since it is derived purely 
on functional grounds. 

3. Functional-type a Posteriori Error Estimates 

In this section, we will first recall the theoretical upper bound for the error in the energy norm, 
and how to minimize this upper bound in order to get a sharp error estimate. Thereafter, we 
will discuss a quality criterion. We will comment on the realization in the isogeometric context 
in Section 4. 

3.1. Guaranteed Upper Bound for the Error. The starting point for the proposed method 
is the following main result, which gives an upper bound for the error in the energy norm. It 
can be found, e.g., in [30, 31, 32]. 



Theorem 3.1. Let Cq be the constant in the FriedricWs type inequality \\v\\ < Cq\\Vv\\a, Vi; E 
Vq. Let u be the exact solution of the problem (4), and let Uh G Vh be an approximate solution. 
Then, the following estimate holds: 

(8) \\Vu - Vu h \\ A < \\AVu h -y\\ A + C n \\ divy + /||, 

where y is an arbitrary vector-valued function in i7(fi,div), and the norms are as defined in (3). 

The constant Cq depends only on the domain ft and the coefficient matrix A (but not on the 
underlying mesh), see, e.g., [25, 32]. Note that Cq can be computed either numerically or, if one 
can find a domain Q\j D fi, where ft\j is a square domain with side-length £, then Cq < C2— -j=., 
where d is the dimension and C2 is the constant in (2). 

Note that, if we choose y via the (unknown) exact solution y = AVu, both sides of (8) coincide. 
Hence, the estimate is sharp in the sense that, for any fixed u^ we can find a function y such 
that the upper bound is as close to the exact error as desired. The estimate given in Theorem 3.1 
is a guaranteed and fully computable upper bound for any conforming approximation Uh G V g . 
In order to obtain a sharp estimate, however, one has to find a function y which minimizes the 
right-hand-side of (8). For minimizing the estimate (8) numerically, we first rewrite the estimate 
in the following form 

(9) \\S7u-Vu h \\ 2 A <(l + P)\\AVu h -y\\ 2 A + (l + ^)C^\\dwy + f\\ 2 =: M%{y,p), 

where /3 > is a free parameter [25, 32]. Note that the upper bound in (9) holds true for any 
fixed y E i?(fi, div) and f3 > 0. Hereinafter, for simplicity, we will refer to M^(y,{3) as the 
majorant. Introducing 



(10) Bi = \\AVu h - y\\% £? 2 = ||divy + /|| 2 , 



B 1 = \\AVi 
we can briefly write the majorant as 

(11) Ml(y,/3) = a 1 B 1 + a 2 B2. 

The efficiency index, defined by 

\\Vu-Vu h \\A 

indicates how close the calculated majorant is to the exact error. The closer I e ^ is to 1, the 
better the estimate. 

It is possible to obtain good estimates by constructing functions y by some post-processing 
of the discrete solution u^, see [25, 32] and the references therein. Consider, for example, that 
Uh E Vh has been computed with NURBS basis functions of degree p > 2. Then, since Uh E C p ~ 1 
(without any repeated internal knot), we have Vuh E (C p ~ 2 ) 2 D i7(fi,div). Choosing y = Vuh 
will thus result in 

(13) ||V«-V« h m<C n ||A« fc + /||. 

To show the efficiency of the estimator (13), we present an illustrative numerical example. This 
example, referred to as Example 1 in the remainder, is chosen due to a smoothly varying function 
(without any large gradients) in both directions. 

Example 1 (Sinus Function on the Unit Square): In this numerical example, the compu- 
tational domain is the unit square ft = (0, l) 2 and Uh is piecewise quadratic in both 
directions, i.e., p = q = 2. The coefficient matrix is constant, A = /, and the exact 
solution is given by 

u — sin(67rx) sin(37ry). 

The right-hand-side / and the (homogeneous) boundary conditions gr> are determined 
by the prescribed exact solution u. 

Once we have calculated t]q := \\Auh + /||q for each cell Q of the mesh, we can compare the 
local errors and choose a criterion for selecting cells which will be marked for further refinement. 
Typically, one chooses a threshold O and marks all cells Q for refinement, where the local error 



is above this threshold. There are several possibilities for determining O, e.g., the bulk-criterion 
proposed in [16]. For simplicity, we choose a percentage i/j and mark a cell Q for refinement, if 

(14) t]q > 6, where 6 = (100 — ^-percentile of {t]q}q. 

The a-percentile of a set A = {ai, . . . ,a^} denotes the value a below which a percent of all 
values ai fall. For example, if we choose ip = 20% in (14), then 6 is chosen such that uq > 6 
holds for 20% of all cells Q. 
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Figure 1. Cells marked by exact error with ^ = 20% in Example 1. 
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Figure 2. Cells marked by error estimator with -0 = 20% in Example 1, y h = Vu^. 
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Figure 3. Convergence of exact error and the majorant (13) for Example 1. 

In Figure 1, we present the cells marked for refinement by the exact error. The cells marked 
for refinement by the majorant given in (13) are presented in Figure 2. We see that starting 
from the mesh 32 x 32, the majorant is able to nicely capture the refinement pattern of exact 
error. However, from a closer look at the convergence of the exact error and the majorant, see 



Figure 3, we find that though such an estimate is a guaranteed upper bound and very cheap to 
compute, it over-estimates the exact error, and its convergence is slower than the exact error 
(due to different operators acting on Uh on both sides). 

In order to obtain a sharp estimate, therefore, it requires to find y E i?(Jl,div) and /? > as 
solutions to the global minimization problem 

(15) min Ml(y,P). 

The technique for finding such minimizing parameters y and f3 will be discussed in Sections 3.2 
and 4.2. Before proceeding further, we give the following Lemma 3.3, which can be found in [32, 
Prop. 3.10]. It provides an analytical result on the sharpness of the bound M^(y,(3). For later 
reference, we also sketch the proof. 

Definition 3.2. A sequence of finite- dimensional subspaces {Yj} ( ^ =1 of a Banach-space Y is 
called limit dense in Y, if for any e > and any v G Y , there exists an index j e , such that 
mf PkeYk \\Pk ~ v\\ Y < £ for all k > j £ . 

Lemma 3.3. Let the spaces {Yj}°? 1 be limit dense in i7(fi,div). Then 

lim inf Ml( yj1 (3) = \\Vu- Vu h \\ 2 A . 

j^oo yjeYj,p>0 

Proof. Recall that the if(fi, div)-norm || • ||di v is defined by |M|^ iv = \\v\\ 2 + || divu|| 2 . Let e > 
be arbitrarily small, but fixed. Let j £ be the index such that, for all k > j £ , there exists a 
Pk G Y k with ||AV^ — PfeHdiv < e - Then, 

(16) inf M*( yj ,P)< M^(p k ,e) = (l + e)\\AVu h -p k \\ 2 A + (l + ±)C{s\\f + div Pk \\ 2 . 

yjeYj,(3>0 

Since ||Ai?||^ = \\v\\a, we can write 

\\AVu h -p k \\ A < \\AVu h - AVu\\ A + \\AVu - p k \\ A 
= \\Vu h -Vu\\ A + \\AVu - p k \\ A . 

The norm || • ||^ is equivalent to the L 2 -norm, so there exists a constant c Al such that the second 
term in the right-hand side can be bounded by 

\\AX7u ~Pk\\A < c A \\AVu ~Pk\\ < c A \\AVu - p k \\&iv < ca£- 
Hence, we obtain the following estimate for the first term in (16): 

(17) \\AVu h - p k \\ A <\\Vu- Vu h \\ A + 0{e). 
Since f — — div AVu, we can bound the second term in (16) as follows: 

(18) \\divp k + f\\ = \\ div p k - div AVu\\<\\p k - AVu\\ div < e. 
With (17) and (18), we can rewrite (16) as 

M|(p fe ,e) < (1 + e)(||Vu - Vu fc |fe + 0{e)) + (1 + ±)C£ £ 2 = ||V« - Vu h \\\ + 0{e). 
Hence, the bound M^(p kl s) — >► \\Vu — Vuh\\ 2 A as e — > 0. D 

3.2. Steps Involved in Minimizing M^(y, /3). As mentioned above, we need to find parame- 
ters y and /3 which minimize the major ant. To do this, we apply an interleaved iteration process 
in which we alternately fix one of the variables and minimize with respect to the other. This 
process, which we summarize in the following, has been described, e.g., in [24, 25]. 

Step 1: Minimization with respect to y: Assume that f3 > is given and fixed, either by an 
initial guess or as a result of Step 2 below. We view the majorant M^(y) as a quadratic 
function of y and calculate its Gateaux-derivative M^(y)' with respect to y in direction 
y. Setting M^{y)' — 0, we obtain 

(19) a\ I A~ y • y dx + a2 / divy div ydx = a\ \ Viih • y dx — a2 / / div ydx, 

Jtt Jn Jn Jn 




where a\ = 1 + (3 and a<i = (1 + t?)Cq, as denned in (10). In order to solve (19), we 
choose a finite-dimensional subspace Y^ C i7(fi,div) and search for a solution y^ E 1^. 
Testing in all directions y ^Y^ leads to a linear system of equations which we write as 

(20) L h y h = r h . 

Here, L h and r h are the matrix and the vector induced by the left hand side and the 
right hand side of equation (19), respectively. By solving (20), we obtain the coefficient 
vector y for the discrete function y^ minimizing M^(y) in Y^ C i7(fi,div). Note that 
this process requires non-negligible cost as we need to assemble L h and r h and solve the 
system (20). 
Step 2: Minimization with respect to f3: Assume that y^ is given from Step 1. By direct 
calculation, we see that M^(/3) is minimized with respect to (3 by setting 

(21) /? = 

where B\ and B^ are as defined in (10). Note that the evaluation of B\ and £>2 (and 
thus /?) requires only the evaluation of integrals, and thus involves negligible cost. 

Steps 1 and 2 are repeated iteratively. We will refer to one loop of applying Step 1 and Step 2 
as one interleaved iteration. Once we have computed minimizers y^ and /3, the computation of 
the majorant M^y^,/?) is straight-forward as it requires only the evaluation of the integrals. 

Note that the matrix L h can be written as 

(22) L h = ai L\ + a 2 L 2 h , 

where L\ and L\ correspond to the terms J^A~ 1 y • y dx and J^divy divy dx in (19), respec- 
tively. Since the matrices L\ and L\ in (22) do not change in the interleaved iteration process, 
they need to be assembled only once. Analogously to (22), we can write r h as 

(23) Lh = a id ~ a 2Ll, 

where r} h and r^ correspond to the terms J n Vu^ • y^ dx and J n f divy dx in (19), respectively. 
The terms r\ and r^ also need to be assembled only once since they also do not change in the 
interleaved iteration process. The full matrix L h and vector r h , however, do change in each 
iteration, because of the change in /3 and y^. Based on past numerical studies, see, e.g., [24, 25], 
and the results presented in Section 5, it has been found that for linear problems, one or two 
such interleaved iterations are enough for obtaining a sufficiently accurate result. 

To recapitulate, we summarize the steps for computing the majorant in Algorithm 1. 

Algorithm 1 Computation of the majorant M e 

Require: u hl f,C^Y h 
Ensure: M e 

P := initial guess 

Assemble and store L^, L^, r^, r^ 

while convergence is not achieved or maximum number of interleaved iterations is not reached 

do 

L h :=(l + ml + (l + ji)C 2 n Li 

Solve L h y h = r h for y_ h 
Bi:= \\AVu h -y h \\\ 

B 2 :=\\di vy h + f \\ 2 

end while 



M e (y, P) := J(l + f3)B 1 + (1 + ±)ClB 2 



i3'^a J 



2 dx. 



3.3. Quality Indicator and Local Error Indicator. So far, we have defined the majorant 
and discussed how we minimize (numerically) the majorant over Y^. Another important ques- 
tion, especially in the light of adaptive, local refinement, is whether a calculated majorant does 
correctly capture the error distribution. From the proof of Lemma 3.3, we recall the following 
observation: 

a\B\ —> \\Vu — Vuh\\A an d ^2^2 -» 0, as y^ E iJ(il,div) — > AVu. 

From this, we deduce the following quality indicator. 

Proposition 3.4. The distribution of the exact error is captured correctly, if 

(24) a 1 B 1 > C© a 2 B 2 

with some constant C© > 1. 

This criterion is easy to check, since the terms appearing in (24) are evaluated in the process 
of minimizing M^(y,(3). It was found in the numerical examples presented in Sections 4 and 5 
that a correct distribution of the error is obtained, if C© > 5. 

We define the local error indicator t]q on a cell Q as the restriction of the majorant to the 
cell Q, i.e., by 

(25) rfc(y h ,P) = (1 + /3) J (Vu h - A- l y h ){AVu h - y h ) dx + (1 + £)C£ J (divy fc + /) 
Obviously the majorant can be expressed as a sum of these local contributions 

Q 

Note that Theorem 3.1 does not guarantee that the local version of the error estimator (25) 
provides a sharp bound for the local error. Hence, some local inaccuracies and over-estimations 
may occur. However, we know that the estimate of the upper bound is guaranteed, and converges 
to the exact error globally. For refinement based on t]q, we again use the criteria (14). 

4. Efficiency and Computational Cost of the Proposed Estimator in the 

isogeometric context 

We now discuss the efficiency and the computational cost of the proposed estimator based 
on the global minimization steps presented in Section 3.2. We again consider Example 1 from 
Section 3.1. All the computations for this example and the examples presented in Section 5 
are performed in matlab®, and the linear systems (5) and (20) are solved using the in-built 
direct solver. One can, however, also use efficient iterative solvers, see, e.g., [19, 20] for (5). 
The right-hand-side / and the boundary conditions go are determined by the prescribed exact 
solution u. 

In the tables, we indicate the mesh-size by the number of interior knot spans of the knot 
vectors s and t, respectively. By this, we mean the number of knot spans without counting the 
vanishing knot spans at the beginning and the end of the open knot vectors. For example, if 

s = (0,0,0,0.25,0.5,0.75,1,1,1) 
t = (0,0,0,0,0.5,1,1,1,1), 

then the mesh-size is 4 x 2. 

We study the efficiency of the majorant based on straight forward computational procedure, 
as discussed in Section 4.1, and based on cost-efficient procedure, as discussed in Section 4.2, 
which coarsens the mesh and increases the polynomial degree simultaneously. This alternative 
cost-efficient procedure will then be used in Section 5 for further numerical examples. In all the 
numerical results of Example 1 in this Section, the initial guess for (3 is 0.01. 

We compare the timings for assembling and for solving the linear systems (5) and (20), as well 
as the total time for assembling and solving. In the presented tables, these timings are shown 
in the columns labelled "assembling-time" , "solving-time" , and "sum" , respectively. The label 
"pde" indicates that the column corresponds to solving the partial differential equation (5), i.e., 
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to assembling K_ h and solving (5) for u h . The label "est." indicates that the timings correspond 
to the estimator, i.e, assembling L h and solving (20) for y . 

In order to check the quality criterion discussed in Section 3.3, we present the values of a\B\ 
and CL2B2 and see whether the inequality (24) is fulfilled or not. To indicate the quality of 
the error distribution captured by the major ant, we plot which cells are marked for refinement 
based on the computed error estimate and the criterion (14). We compare this to the refinement 
marking based on the criterion (14) applied to the exact local error instead of the estimate. 

4.1. Straightforward Procedure. Analogously to Vh and Vh, we choose a function space Yh 
on the parameter domain and we define the function space Yh by the push-forward 

Y h = Y h oG~\ 

For our first choice for Y^, we use the same mesh as for Vh, and we choose 



(26) 



Yh — Sp+l^q ® <Sp } g+l, 



where S Piq denotes the space of NURBS functions of degree p and C p_ ^continuity in the first 
coordinate, and degree q and C g_ ^continuity in the second coordinate. 

In Table 1, we present the computed efficiency indices obtained with this choice of Yh, which 
show that sharp upper bounds are obtained as the mesh is refined. The dashed line in Table 1 
indicates that the criterion (24) is fulfilled with C& = 2 starting from the mesh 32 x 32. 

The cells marked by the error estimator are shown in Figure 4, plotted in magenta. When 
comparing these plots, we see that the error distribution is captured accurately starting from 
the mesh 32 x 32. 

The timings presented in Table 2, however, show that the computation of the error estimate 
is costlier (about 4.5 times) than assembling and solving the original problem. This is not 
surprising, since, when N u denotes the number of degrees of freedom (DOF) of Uh, the number 
of DOF of yh, which is vector- valued, is 2N U (asymptotically) and the linear system is solved 
with a direct solver. Clearly, this straightforward approach is not cost-efficient. 



Table 1. 

in (26). 



mesh-size 


htt 


a 1 B 1 


a 2 B 2 


8x8 


3.43 


2.62e+01 


1.17e+02 


16 x 16 


1.92 


6.07e-01 


6.19e-01 


32 x32 


1.41 


2.29e-02 


9.71e-03 


64x64 


1.20 


1.15e-03 


2.33e-04 


128 x 128 


1.10 


6.51e-05 


6.54e-06 


256 x 256 


1.05 


3.87e-06 


1.95e-07 


512 x 512 


1.03 


2.36e-07 


5.94e-09 



Efficiency index and components of the majorant in Example 1, Yh as 
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■H"l 1- 



(a) 16 x 16 



M I II 
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M I M 

(b) 32 x 32 
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• •III 

• MM 

(c) 64 x 64 



• MM 
IMM 

• MM 

(d) 128 x 128 



Figure 4. Cells marked by error estimator with ip = 20% in Example 1, Yh as in (26). 
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mesh-size 


#DOF 


assembling-time 


solving-time 


sum 




Uh 


Vh 


pde 


est. 


u est. •>•> 
pde 


pde est. 


a est, rj 
pde 


pde est. "^" 


8x8 


100 


220 


0.04 


0.17 


4.39 


<0.01 <0.01 


5.16 


0.04 0.17 4.40 


16 x 16 


324 


684 


0.14 


0.59 


4.25 


<0.01 0.01 


5.39 


0.14 0.60 4.26 


32 x32 


1156 


2380 


0.46 


2.17 


4.70 


0.01 0.03 


4.71 


0.47 2.20 4.70 


64x64 


4356 


8844 


1.82 


8.51 


4.68 


0.03 0.20 


6.15 


1.85 8.70 4.70 


128 x 128 


16900 


34060 


7.38 


34.19 


4.63 


0.15 0.87 


5.70 


7.54 35.06 4.65 


256 x 256 


66564 


133644 


33.30 


149.78 


4.50 


0.84 5.66 


6.78 


34.14 155.44 4.55 


512 x 512 


264196 


529420 


191.11 


766.10 


4.01 


3.77 33.92 


9.00 


194.88 800.03 4.11 




Table 2. Number of DOF and timin 


gs in Example 


l,Y h i 


is in (26). 



4.2. Alternative Cost-Efficient Procedure. Recall that the cost of Step 1 of the algorithm 
presented in Section 3.2 depends on the choice of Y^ C if(f2, div). As shown in Lemma 3.3, we 
can make the estimate as sharp as we desire by choosing a suitably large space Y^. However, the 
larger Y^ is chosen, the more costly setting up and solving the system (20) becomes. Clearly, it 
is highly desirable to keep the cost for error estimation below the cost for solving the original 
problem. 

As discussed above, choosing Y^ as in (26) does not result in a cost-efficient method. Apart 
from the fact that y^ is vector-valued while u^ is scalar, another aspect contributes to the high 
cost for the procedure presented in Section 4.1. Recall that, by choosing Y^ as in (26), we have 

y\ E 5p+i, g , 

i.e., the components of y^ are in different spline spaces. Hence, we have to compute different 
basis functions for y\ and y<i (note that this can be a costly procedure for higher polynomial 
degrees). Furthermore, when assembling, for example, the matrix L h , we need to compute 
integrals over products of basis functions of the form 

RiRj dx. 

With Yh as in (26), the product RiRj of basis functions of y\ is different to the product of basis 
functions of 3/2? hence, the integrals have to be evaluated independently for y\ and 7/2 • 
In the light of these observations, we propose the following alternative choice for Y^. 



(27) 



Yh = Sp+l,q+l ® <Sp+l j9 +l. 



We refer to this setting as Case 1 in the remainder of the paper. With this choice, y\ and y<i 
are contained in the same spline spaces. Hence, the basis functions need to be computed only 
once, and any computed function values can be used for both components of y^. 

The computed efficiency indices are presented in Table 3, which show that we obtain even 
better (i.e., sharper) upper bounds for the exact error with Y^ as in (27) than with the choice 
(26). When we compare the plots of the cells marked by the error estimator in Figure 5 to the 
plots in Figure 1, we see that the error distribution is again captured accurately starting from 
the mesh 32 x 32. The dashed line in Table 3 indicates that the criterion (24) is fulfilled with 
C0 = 3 starting from the mesh 32 x 32. 

The timings obtained with this method are presented in Table 4. This approach reduced the 
total time needed for computing the majorant from a factor of about 4.5 to a factor of approx- 
imately 3 compared to the time for assembling and solving the original problem. Nevertheless, 
a factor of 3 in the timings is still not satisfactory. 

In order to further reduce the computational cost, we reduce the number of DOF of y^ by 
coarsening the mesh by a factor K in each dimension. The number of DOF of y^ is thus reduced 
to 2N U /K 2 (asymptotically). The larger K is chosen, the greater the reduction of DOF will be. 
At the same time, if the coarsening is done too aggressively, sharp features might not be detected 
properly on coarse meshes. We counter the reduction in accuracy due to mesh-coarsening by 
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mesh-size 


-feff 


a 1 B 1 


a 2 B 2 


8x8 


2.77 


8.08e+01 


1.24e+01 


16 x 16 


1.71 


5.75e-01 


3.96e-01 


32 x32 


1.32 


2.14e-02 


7.05e-03 


64x64 


1.16 


l.lle-03 


1.78e-04 


128 x 128 


1.08 


6.39e-05 


5.08e-06 


256 x 256 


1.04 


3.83e-06 


1.53e-07 


512 x 512 


1.02 


2.35e-07 


4.69e-09 



Table 3. Efficiency index and components of the majorant in Example 1, Case 1. 



HbH 


i ■■■ i 


H"H 



(a) 16 x 16 



MIM 

tint 

MIM 

(b) 32 x 32 



IMHII 


• 1*11 


IMHII 



(c) 64 x 64 



• •Ml 
IMM 

• •••• 

(d) 128 x 128 



Figure 5. Cells marked by error estimator with if; = 20% in Example 1, Case 1. 



mesh-size 


#DOF 


assembling-time 


solving-time 


sum 




Uh 


Vh 


pde 


est. 


u est. •>•> 
pde 


pde 


est. 


u est. •>•> 
pde 


pde est. "^" 


8x8 


100 


242 


0.04 


0.11 


2.78 


<0.01 


<0.01 


1.51 


0.04 0.11 2.76 


16 x 16 


324 


722 


0.12 


0.34 


2.86 


<0.01 


0.01 


5.33 


0.12 0.35 2.90 


32 x32 


1156 


2450 


0.46 


1.35 


2.94 


0.01 


0.05 


7.69 


0.47 1.40 3.01 


64x64 


4356 


8978 


1.77 


5.30 


2.99 


0.03 


0.27 


8.02 


1.80 5.57 3.09 


128 x 128 


16900 


34322 


7.39 


21.89 


2.96 


0.16 


1.45 


9.26 


7.55 23.34 3.09 


256 x 256 


66564 


134162 


33.00 


94.69 


2.87 


0.84 


8.83 


10.54 


33.84 103.52 3.06 


512 x 512 


264196 


530450 


191.59 


498.20 


2.60 


3.83 


61.45 


16.06 


195.42 559.65 2.86 



Table 4. Number of DOF and timings in Example 1, Case 1. 



increasing the polynomial degree of y^ by some positive integer &;, i.e., by choosing 

(28) Yh = S p +k,q+k ® <Sp+k,q+k- 

Note that one could also choose different factors K\ and K 2 and different degree increases k\ 
and k 2 for the first and second component, respectively, if desired. 

Remark 4.1. With these choices ofY^, we take advantage of the following specific property of 
univariate NURBS basis functions. For C p ~ 1 regularity, increasing the polynomial degree by k 
only adds a total of k additional basis functions. In other words, the global smoothness can be 
increased at the cost of only a few additional DOF. Coarsening the mesh by a factor K , however, 
will also reduce the number of DOF by the same factor K (asymptotically) . 

Note that Case 1 discussed above fits into this framework, since Case 1 corresponds to the 
choice K — k — 1. For the next setting, we apply moderate mesh-coarsening by choosing 

K = k = 2 (i.e., Y h = S p+2 , q +2 ® S p+2 ^ q+2 ). 

This setting will be referred to as Case 2 in the remainder of the paper. 

The computed efficiency indices along with the magnitudes of the terms a\B\ and a 2 B 2 for 
Case 2 are presented in Table 5, and the marked cells are plotted in Figure 6. Both indicate 
that a good upper bound of the error and the correct error distribution are computed on fine 
meshes. On coarse meshes, however, the efficiency index is larger than in Case 1, which is due 
to the mesh-coarsening. 
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The timings presented in Table 6 show that, even though Case 2 is faster than Case 1, this 
approach still costs roughly as much as solving the original problem. This is due to the costlier 
evaluation of the higher degree basis functions, as well as the increased support and overlap of 
the basis functions, which results in more non-zero entries in L h than in K h . 



mesh-size 


left 


a\B\ a 2 B 2 


8x8 


14.19 


1.59e+03 8.53e+02 


16 x 16 


8.49 


1.97e+01 4.32e+00 


32 x32 


1.82 


3.05e-02 2.41e-02 


64 x 64 


1.16 


1.12e-03 1.76e-04 


128 x 128 


1.04 


6.14e-05 2.24e-06 


256 x 256 


1.01 


3.72e-06 3.32e-08 


512 x 512 


1.00 


2.31e-07 5.13e-10 



Table 5. Efficiency index and components of the majorant in Example 1, Case 2. 




M*M 



***** 



If til 
11*14 
***** 



MMI 
Mill 

MMI 



(a) 16 x 16 (b) 32 x 32 (c) 64 x 64 (d) 128 x 128 

Figure 6. Cells marked by error estimator with ^ = 20% in Example 1, Case 2. 



mesh-size 


#DOF 


assembling-time 


solving-time 




sum 




Uh 


Vh 


pde 


a «est. n 
ebT; - pde 


pde est. "f|" 


pde 


a «est.» 
ebt - pde 


8x8 


100 


128 


0.03 


0.05 1.39 


<0.01 <0.01 1.16 


0.04 


0.05 1.39 


16 x 16 


324 


288 


0.14 


0.18 1.29 


<0.01 <0.01 0.92 


0.14 


0.18 1.28 


32 x 32 


1156 


800 


0.54 


0.59 1.10 


0.01 0.02 2.32 


0.55 


0.61 1.11 


64 x 64 


4356 


2592 


1.91 


2.33 1.22 


0.04 0.08 2.09 


1.95 


2.40 1.23 


128 x 128 


16900 


9248 


7.46 


9.54 1.28 


0.19 0.51 2.75 


7.64 


10.05 1.32 


256 x 256 


66564 


34848 


33.93 


39.02 1.15 


0.90 2.59 2.88 


34.82 


41.60 1.19 


512 x 512 


264196 


135200 


196.23 


177.98 0.91 


4.08 15.91 3.90 


200.31 


193.89 0.97 




Table 6. N 


umber < 


Df DOF and tii 


mings in Example 1, 


Case 2. 





In order to further improve the timings, we coarsen the mesh more aggressively by a factor 
of 4 and, at the same time, increase the polynomial degree of yh by 4 compared to Uh, i.e., 

K = k = 4 (i.e., Y h = 5 p+4 , g +4 ® Sp+4,g+4). 

We refer to this setting as Case 3 in the remainder of the paper. 

This aggressive coarsening notably affects the efficiency index on coarse meshes, see Table 7. 
On fine meshes, however, the efficiency indices are close to 1 in all presented cases. The number 
of DOF of yh in Case 3 is only N u /8 (asymptotically). The timings presented in Table 8 show 
that this setting results in a method which can be performed significantly faster (at almost half 
of the cost) than solving the original problem. The more aggressive reduction of DOF outweighs 
the additional costs mentioned above, even though the polynomial degree is now increased by 4. 

Remark 4.2. In all cases for Example 1, criterion (24) is fulfilled with C^ = 5 on meshes of 
size 64 x 64 and finer. This is indicated by the dashed lines in Tables 5 and 7 for Cases 2 and 3, 
respectively, and in Figures 6 and 7. As we have discussed above, the error distribution might 
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already be captured with a smaller constant C^ in criterion (24), but Cases 2 and 3 show that 
this is not guaranteed. Example 1, and also the examples discussed in Section 5, indicate that 
Cq = 5 is a good choice for checking criterion (24) numerically, even though this choice may be 
conservative in some cases. 



mesh-size 


left 


a\B\ a 2 B 2 


8x8 


11.28 


5.38e+02 1.01e+03 


16 x 16 


36.43 


2.83e+02 1.60e+02 


32 x32 


12.63 


2.04e+00 5.81e-01 


64 x 64 


1.17 


1.13e-03 1.88e-04 


128 x 128 


1.01 


5.98e-05 3.79e-07 


256 x 256 


1.00 


3.70e-06 1.24e-09 


512 x 512 


1.00 


2.31e-07 5.32e-12 



Table 7. Efficiency index and components of the majorant in Example 1, Case 3. 



■ 1 






■ ■ 
■ 


■ ■ 
■ 


l| 


| j 




il 

1 


li 

1 


rf 


li 




■ I 

■ 

■ ■ 


|i 

■ 
■ 1 



HIM 
• If •• 
MtM 



• •••I 
MM! 

• •••I 



(a) 16 x 16 (b) 32 x 32 (c) 64 x 64 (d) 128 x 128 

Figure 7. Cells marked by error estimator with ip = 20% in Example 1, Case 3. 



mesh-size 


#DOF 


assembling- 


time 


solving-time 




sum 




Uh 


Vh 


pde 


est. 


(test, i) 
pde 


pde est. 


uest. •>•) 
pde 


pde 


„j. «est. » 
ebl - pde 


8x8 


100 


128 


0.04 


0.03 


0.76 


<0.01 <0.01 


1.09 


0.04 


0.03 0.76 


16 x 16 


324 


200 


0.14 


0.10 


0.69 


<0.01 <0.01 


0.61 


0.14 


0.10 0.69 


32 x32 


1156 


392 


0.54 


0.31 


0.57 


0.01 <0.01 


0.34 


0.55 


0.31 0.57 


64 x 64 


4356 


968 


1.90 


1.19 


0.63 


0.04 0.01 


0.26 


1.94 


1.20 0.62 


128 x 128 


16900 


2888 


7.49 


4.86 


0.65 


0.16 0.14 


0.84 


7.66 


4.99 0.65 


256 x 256 


66564 


9800 


33.90 


20.15 


0.59 


0.91 0.82 


0.91 


34.81 


20.98 0.60 


512 x 512 


264196 


35912 


194.25 


84.70 


0.44 


4.10 5.45 


1.33 


198.35 


90.15 0.45 



Table 8. Number of DOF and timings in Example 1, Case 3. 



Remark 4.3. The observations discussed above illustrate that one has to balance the sharpness 
of the majorant on the one hand, and the required computational effort on the one hand. Note 
that in typical practical applications, the exact solution (and thus the sharpness of the majorant) 
is not known. Therefore, to address the balancing between sharpness and required computational 
effort, we propose the following strategy. If the mesh is coarse and the total computational cost 
for the error estimate is still moderate, we apply no (or only moderate) coarsening. When the 
original mesh is already fine, that is, if the problem size is large, we coarsen the mesh more 
aggressively and thereby profit from the fast computation of the estimate. 

We now comment on the interleaved iterations. The results in the Tables 1, 3, 5, and 7 were 
obtained by applying only two interleaved iterations, as described in Section 3.2. As mentioned 
there, a sufficiently accurate result can be obtained already after the first such iteration. To 
illustrate this, we present the efficiency indices for Case 3 in Table 9, which were obtained after 
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one, two, and four interleaved iterations, respectively. The efficiency index does vary notably 
on the coarser meshes, but since all of these values greatly overestimate the exact error, they 
do not correctly capture the error distribution. On meshes, where the criterion (24) is fulfilled, 
and thus the error distribution is correctly recovered, the differences due to more interleaved 
iterations are insignificant. 



mesh-size 


interleaved iterations 




1 2 


4 


8x8 


11.84 11.28 


11.25 


16 x 16 


80.31 36.43 


33.78 


32 x32 


17.36 12.63 


10.11 


64 x 64 


1.20 1.17 


1.17 


128 x 128 


1.01 1.01 


1.01 


256 x 256 


1.00 1.00 


1.00 


512 x 512 


1.00 1.00 


1.00 


of I e flf for different numbi 


3rs of int 



ample 1, Case 3. 



5. Numerical Examples 

In this section, we present further numerical examples which illustrate the potential of the 
proposed a posteriori error estimator. We will discuss the following three settings that were also 
discussed in Section 4. Case 1: K = k = 1, Case 2: K = k = 2, and Case 3: K = k = 4. As in 
Example 1, the initial guess for f3 is 0.01. 

The parameter domain in all presented examples is the unit square Cl = (0, l) 2 , as discussed 
in Section 2.2. The mesh-sizes in the two coordinate directions, which will be presented in the 
tables, are determined by the respective initial meshes, which in turn, are determined by the 
geometry mappings. The figures plotted in black represent the computations based on the exact 
error, and the figures plotted in magenta represent the computations based on the majorant. 





(a) Example 2. a (a = 20). 



(b) Example 2.b (a = 50). 



Figure 8. Exact solutions u on. Q. 

Example 2: Quarter Annulus. The domain Q for our second example represents a quarter 
annulus. In polar coordinates, fi is defined by (r, (j)) E (1, 2) x (0, ^). Note that the circular parts 
of the domain boundary are represented exactly by the NURBS geometry mapping of degree 2, 
i.e., we have p = q = 2. We set A = I, and we prescribe the exact solution 

u={r- l)(r - 2)<t>{(t> - z) e -*(rcos<p-i)\ 

We test our method with two values of a, namely, 

Example 2. a: a = 20, Example 2.b: a = 50. 

In both examples, this function has zero Dirichlet boundary values and a peak at x = 1, the 
sharpness of which is determined by the value of a. The exact solutions are depicted in Figure 8. 
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This example is chosen because of sharp peaks (large gradients) and curved boundary (requiring 
NURBS mapping for exact representation). 



mesh-size I e Q a\B\ CL2B2 





Case 1 


16x8 


1.83 


9.98e-04 3.59e-04 


32 x 16 


1.29 


2.08e-05 6.51e-06 


64x32 


1.13 


1.04e-06 1.44e-07 


128 x 64 


1.07 


5.95e-08 4.00e-09 


256 x 128 


1.03 


3.58e-09 1.20e-10 


512 x 256 


1.02 


2.20e-10 3.67e-12 



Case 2 



16 x 8 


13.99 


4.44e-02 3.51e-02 


32 x 16 


4.17 


2.00e-04 8.43e-05 


64x32 


1.31 


1.20e-06 3.66e-07 


128 x 64 


1.06 


5.91e-08 3.36e-09 


256 x 128 


1.01 


3.51e-09 4.60e-ll 


512 x 256 


1.00 


2.17e-10 6.96e-13 





Case 3 


16 x 8 


24.87 


1.09e-01 1.42e-01 


32 x 16 


56.02 


2.92e-02 2.22e-02 


64x32 


10.42 


7.81e-05 2.16e-05 


128 x 64 


1.11 


6.21e-08 6.61e-09 


256 x 128 


1.00 


3.49e-09 1.02e-ll 


512 x 256 


1.00 


2.17e-10 3.27e-14 



Table 10. Efficiency index and components of the majorant in Example 2. a (a = 20). 




(a) Exact, 
mesh 32 x 16. 




(e) Exact, 
mesh 64 x 32. 




(i) Exact, 
mesh 128x64. 




(b) Case 1, 
mesh 32 x 16. 




(f) Case 1, 
mesh 64 x 32. 




(j) Case 1, 
mesh 128x64. 




(c) Case 2, 
mesh 32 x 16. 




(g) Case 2, 
mesh 64 x 32. 



^ 



(k) Case 2, 
mesh 128x64. 




(d) Case 3, 
mesh 32 x 16. 




(h) Case 3, 
mesh 64 x 32. 




(1) Case 3, 
mesh 128x64. 



Figure 9. Marked cells with ip = 20% in Example 2. a (a = 20). 

In Tables 10 and 11, the efficiency index I e Q and the magnitudes of a\B\ and OL2B2 are 
presented for both examples. The dashed lines indicate the mesh-size after which criterion (24) 
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mesh-size I-a a\B\ CL2B2 





Case 1 




16x8 


3.02 


2.94e-02 


1.78e-02 


32 x 16 


1.92 


3.57e-04 


1.83e-04 


64x32 


1.34 


9.15e-06 


3.22e-06 


128 x 64 


1.16 


4.67e-07 


7.56e-08 


256 x 128 


1.08 


2.67e-08 


2.12e-09 


512 x 256 


1.04 


1.60e-09 


6.32e-ll 





Case 2 




16 x 8 


13.84 


3.45e-01 


6.49e-01 


32 x 16 


16.76 


2.58e-02 


1.53e-02 


64x32 


3.16 


4.10e-05 


2.80e-05 


128 x 64 


1.25 


5.04e-07 


1.24e-07 


256 x 128 


1.05 


2.61e-08 


1.33e-09 


512 x 256 


1.01 


1.56e-09 


1.89e-ll 



Case 3 



16 x 8 


17.20 


4.24e-01 l.lle+00 


32 x 16 


76.95 


3.24e-01 5.41e-01 


64x32 


83.72 


3.02e-02 1.83e-02 


128 x 64 


4.19 


4.64e-06 2.44e-06 


256 x 128 


1.04 


2.59e-08 1.02e-09 


512 x 256 


1.00 


1.55e-09 2.22e-12 



Table 11. Efficiency index and components of the majorant in Example 2.b (a = 50). 




^^ 




(a) Exact, 
mesh 64 x 32. 



(b) Case 1, 
mesh 64 x 32. 





(e) Exact, 
mesh 128x64. 



(f) Case 1, 
mesh 128x64. 



(c) Case 2, 
mesh 64 x 32. 



^ 



(g) Case 2, 
mesh 128x64. 



(d) Case 3, 
mesh 64 x 32. 




(h) Case 3, 
mesh 128x64. 




^^% 



(i) Exact, 
mesh 256 x 
128. 



(j) Case 1, 
mesh 256 x 
128. 



(k) Case 2, 
mesh 256 x 
128. 



(1) Case 3, 
mesh 256 x 
128. 



Figure 10. Marked cells with ip = 20% in Example 2.b (a = 50). 

with C© = 5 is fulfilled. The distribution of the marked cells is depicted in Figures 9 and 10. 
As before, we observe that the error distribution is represented correctly if the criterion (24) is 
fulfilled with C© = 5. 

When comparing Tables 10 and 11, as well as Figures 9 and 10, we notice the following. 
The more aggressive the mesh coarsening, and the sharper the peak, the more refinements are 
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needed before criterion (24) is fulfilled and the error distribution is captured correctly. Also, the 
obtained efficiency index, and thus, the quality of the error bound, is better when only moderate 
coarsening is applied. 

Since the timings in Example 2. a and Example 2.b show the same behaviour as in Example 1, 
both regarding assembling-time and solving-time, we omit the presentation of these numbers. 
Clearly, Case 3 outperforms Cases 1 and 2 in terms of cost-efficiency. 

Example 3: L-shaped Domain. In our third example, we consider the Laplace equation 

(29) Au = 

with Dirichlet boundary conditions on the L-shaped domain £1 = (— 1, 1) 2 \[0, l] 2 . As before, 
p = q = 2. The function 

2 

u(r, (ft) = r 3 sin((20 - 7r)/3) 

solves (29) and is used to prescribe Dirichlet boundary conditions. The solution has a singularity 
at the re-entrant corner at (0,0). This is a classical example for a posteriori error estimation. 



mesh-size J e ff a\B\ CL2B2 



Case 1 



32 x 16 


2.73 


2.94e-03 2.87e-04 


64 x 32 


2.54 


1.09e-03 1.01e-04 


128 x 64 


2.47 


4.39e-04 3.75e-05 


256 x 128 


2.49 


1.87e-04 1.50e-05 





Case 2 


32 x 16 


3.82 


5.68e-03 6.44e-04 


64 x 32 


3.39 


1.91e-03 2.09e-04 


128 x 64 


3.03 


6.47e-04 6.74e-05 


256 x 128 


2.82 


2.37e-04 2.29e-05 





Case 3 


32 x 16 


5.72 


1.28e-02 1.40e-03 


64 x 32 


5.42 


4.87e-03 5.20e-04 


128 x 64 


4.53 


1.45e-03 1.51e-04 


256 x 128 


3.78 


4.22e-04 4.46e-05 



Table 12. Efficiency index and components of the majorant in Example 3. 



HL 



ir 



(a) Exact error. 



(b) Case 1. 



(c) Case 2. 



v 



(d) Case 3. 



FIGURE 11. Marked cells with i/j = 10% in Example 3 on the mesh 128 x 64. 

The magnitudes of the components a\B\ and <22£>2, which are presented in Table 12, show that 
criterion (24) with C§ = 5 is fulfilled on all the considered meshes. As an example, the marked 
cells based on the exact solution and on the error estimator are presented in Figure 11. They 
are plotted for a mesh of size 128 x 64 and show that the singularity near the corner is captured 
by the error estimate. Since u ^ i7 2 (fi), increasing the degree of approximation functions 
(NURBS) contaminate the majorant quality. Therefore, the efficiency indices in this example 
are larger than in the previous examples. The associated cost exhibits a similar behaviour as in 
the previous examples, and is thus not reported here. 
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Example 4: Advection-Dominated Advection-Diffusion-Equation. In our fourth exam- 
ple, we consider the advection-diffusion equation with Dirichlet boundary conditions on the unit 
square £1 = (0, l) 2 , with p = q = 2, i.e., 

—kAu + b • Vu = in fi, 

u — go on <9fi, 



where 



k= 10" 



b = (cos ? , sin ? 



7T\T 

3/ ' 



to = 



1, ifi/ = 

0, else 



With these parameters, the equation is ad vect ion-dominated. The motivation for this example 
is due to the sharp boundary layers. 

We use the standard streamline upwind Petrov-Galerkin (SUPG) scheme for the stabilization. 
The stabilization parameter r is set to r(Q) = /i&(Q)/2|&|, where h^Q) is the diameter of the 
cell Q in direction of the flow 6, and |6| is the magnitude of the vector b. For advection-diffusion 
problems, we have to adapt the majorant. Since the principle method is the same, we refer the 
reader to [ , Section 4.3.1] for a detailed discussion. In this special case, where A — nl with 
k <C |6|, and with constant velocity vector 6, the majorant M { 
problem is given by 



2 adv f° r ^ ne advection-diffusion 



Mi 



,adv 



(l + /3)\\AVu h -y 



W + d 



4)Cg|| div j/ + / - 6 • Vu fc || 5 



P> 



The strong advection and the discontinuous boundary conditions result in sharp layers. 
Figure 12(a), the expected positions of the layers are indicated by dashed lines. 



In 



mesh-size a\B\ (12B2 



Case 1 



16 x 16 


1.98e-07 3.18e-10 


64 x 64 


6.45e-07 1.15e-09 


256 x 256 


2.28e-06 4.33e-09 





Case 2 


16 x 16 


1.83e-06 9.66e-10 


64 x 64 


6.50e-06 3.65e-09 


256 x 256 


1.86e-05 1.24e-08 



Case 3 



16 x 16 


3.24e-06 1.29e-09 


64 x 64 


2.07e-05 6.52e-09 


256 x 256 


6.86e-05 2.38e-08 



Table 13. Comparison of terms a\B\ and CL2B2 in Example 4. 






(a) Expected posi- 
tions of sharp layers. 



(b) Marked cells with 
ip = 20%, mesh-size 
64 x 64. 



(c) Marked cells with 
ip = 10%, mesh-size 
256 x 256. 



Figure 12. Expected layers and marked cells in Example 4, Case 3. 

The magnitudes of a\B\ and CL2B2 presented in Table 13 indicate that the criterion (24) with 
C0 = 5 is fulfilled on all the considered meshes. The distribution of the marked cells presented 
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mesh-size 



#DOF 

Uh Vh 



assembling- 
pde est. 



time 

u est. •>•> 
pde 



solving-time 
pde est. "gk" 



pde 



sum 
est. 



i est. r> 
pde 



Case 1 



16 x 16 


324 722 


0.25 0.39 1.56 


<0.01 


0.01 


6.38 


0.25 0.40 1.59 


64x64 


4356 8978 


3.25 5.32 1.63 


0.03 


0.26 


8.64 


3.28 5.58 1.70 


256 x 256 


66564 134162 


51.22 94.15 1.84 


0.85 


8.84 


10.35 


52.07 102.99 1.98 



Case 2 



16 x 16 


324 


288 


0.21 0.14 0.67 


<0.01 <0.01 


0.50 


0.21 


0.14 0.67 


64 x 64 


4356 


2592 


3.26 2.10 0.64 


0.03 0.06 


2.01 


3.29 


2.16 0.66 


256 x 256 


66564 


34848 


50.83 35.58 0.70 


0.85 2.30 


2.70 


51.68 


37.87 0.73 



Case 3 



16 x 16 


324 


200 


0.26 0.10 0.39 


<0.01 <0.01 


0.58 


0.26 


0.10 0.40 


64 x 64 


4356 


968 


3.41 1.21 0.35 


0.04 0.01 


0.26 


3.44 


1.22 0.35 


256 x 256 


66564 


9800 


52.40 19.83 0.38 


1.02 0.91 


0.89 


53.42 


20.74 0.39 



Table 14. Timings in Example 4. 



in Figures 12(b) and 12(c) provides the visual indication that the expected layers are detected 
by the error estimate. 

In Table 14, the timings are presented. Note that, unlike the previous examples, assembling 
and solving the system for the estimator is faster than for the original problem not only in Case 3, 
but also in Case 2. This is due to the SUPG stabilization which is costlier than computing the 
additional term b • Vu^ in the major ant M^ adv . 

6. Conclusion 

We have proposed a method for cost-efficient computation of guaranteed and sharp a posteriori 
error estimates in IGA. This method relies only on the use of NURBS basis functions, without the 
need for constructing complicated basis functions of H(£l, div). Two properties of NURBS basis 
functions are exploited. Firstly, the basis functions are, in general, automatically in i7(fi,div) 
due to their high smoothness. Without this property, we could not use NURBS as basis functions 
for the minimizing function y^. Secondly, increasing the polynomial degree of NURBS basis 
functions adds only few basis functions. This fact is necessary for keeping the computational 
cost of the majorant as low as possible (see Remark 4.1). 

We have discussed different settings which allow the user to balance the sharpness of the 
bound and accurate error distribution on the one hand, and the required computational cost of 
the error estimator on the other hand (see Remark 4.3). For the presented settings, we have 
derived a quality criterion, which is easy to check and which indicates whether the computed 
estimate is sharp or not (see Remark 4.2). 

In this paper, we have only considered tensor-product NURBS discretizations. While the 
extension of this method to locally refined isogeometric discretizations and also to three dimen- 
sions is, in theory, straightforward, the actual performance and efficiency of the error estimator 
on such methods and meshes is the subject of further studies. 
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